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^N| ■ Phase transition in fullerenes Cgo and C240 are investigated by means of constant- 

temperature molecular dynamics simulations. In the phase transition region, the 
00 ' assembly (and fragmentation) of the Cgo cage from (and to) the gaseous state is 

demonstrated via the dynamical coexistence of two phases. In this critical region, 
^ . the fullerene system is seen to continuously oscillate between the carbon cage (the 

^ ' solid phase) and the state of carbon dimers and short chains (the gas phase) . These 

oscillations correspond to consecutive disintegration and formation of the fullerene. 
c^ ■ 

(713 ' Furthermore, the temperature-dependent heat capacity of the fullerene features a 

O 
'^ , prominent peak, signifying the finite system analogue of a first-order phase transi- 

>^ 

,^ \ tion. The simulations were conducted for 500 ns using a topologically-constrained 



Dh. 



pairwise forcefield which was developed for this work. Results of the simulations 



^ . were supplemented by a statistical mechanics analysis to account for entropy and 

ff^ , pressure corrections, corresponding to experimental conditions. These corrections 

t:J" ' lead to a phase transition temperature of 3800-4200 K for pressure 10-100 kPa, in 

f— N ■ good agreement with available experimental values. 
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I. INTRODUCTION 

Thermal fragmentation and assembly of fullerenes can be viewed as reverse processes 
which have clear features of phase transition. The former leads to the disintegration of the 
solid-like hollow cage to a gas-like state of dimers (C2 units); while the latter recreates the 
fullerene cage from the hot carbon gas. This work investigates the above idea by means 
of isothermal molecular dynamics (MD) simulations of fullerenes Cgo and C24o- We report 
the occurrence of fullerene sublimation above a critical temperature — the phase transition 
temperature — beginning with the loss of a C2 unit. This rapidly leads to the breaking of the 
fullerene cage and its eventual decomposition into the gaseous phase. In the region of the 
phase transition temperature, the system oscillates between two distinct phases: the solid- 
like cage and the gas-like state of the carbon dimers. Such oscillation corresponds to the 
consecutive back-and-forth fragmentation and re-assembly of the fullerene which can be seen 
clearly in the bimodal distribution of the total energy over time (at a single temperature). 
This coexistence behavior, as well as the prominent peak of the temperature-dependent heat 
capacity are signatures of first-order phase transition in__finite systems (l| . 
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Whether the fragmentation is induced by pyrolysis [21, laser-irradiation Sil or collisions 

nnn 

with charged/neutral particles [4, 15|, la] , fullerenes are known to disintegrate via sequential 
loss of C2 units, through asymmetric fission or multifragmentation 7|, |8|. Fullerene sta- 
bility and its fragmentation mechanism have been investigated extensively by a variety of 
computational methods, notably tight-binding molecular dynamics (TBMD). Initial TBMD 
simulations of the Ceo fragmentation were performed by Wang et al. in which the Ceo was 
found to be stable against spontaneous disintegration for temperatures up to 5000 K J9l]. 
This was followed by similar studies by Zhang et al. for fullerenes ranging from C20 to 
C90 where, for small fullerenes (n<58) it was discovered that the fragmentation tempera- 
ture increased linearly. However, for larger fullerenes (n=60 and n>70), the temperature 
stabilised around 5500 K |lOl . Ill| . Both studies employed the TB-parametrisation of Xu 



et al. [I2I], as are the work conducted by Laszlo [13|] and Openov and Podlivaev [ij] for 
the canonical and microcanonical ensembles respectively. Kim et. al reported structural 
changes in the Ceo and C70 in the range of 3000 and 4000 K — with the onset of bond break- 
ing around 5000 K [151]. Kim and Tomanek conducted fragmentation simulations of the C20, 
Ceo and C240 and found several different phases of the fullerene melting process, including a 



liquid-like pretzel phase [l6|. While Xu and Scuseria conducted photofragmentation simula- 
tions of the Ceo and observed sequential loss of C2 units with cage fragmentation occurring 
for T>5600 K. Horvath and Beu have also conducted radiation-induced fragmentation of 
fuUerenes and reported multifragmentation to be the main disintegration channel at high 
excitation energies. However, between excitation energies of 100-120 eV, they reported the 
occurrence of a phase transition where they had defined phase transition to be a steep drop 
in the average fragment size with temperature (in their work, the average fragment size 
dropped from 60 to 5 between 100-120 eV) [17]. Semiempirical bond-order methods are 
also used in studying fuUerenes. The many-body Tersoff potential [18| has been applied to 
a large number of carbon systems, and has been used by Marcos, et al. to investigate the 
thermal stability of fuUerenes [19|. The Reactive Bond Order Potential (REBO, or Bren- 
ner's potential) [20!] — itself built upon the Tersoff potential — is also typically employed and 
has been used to study the formation of carbon-cage structures [2l| and it collision-induced 
fragmentation 22[. 

FuUerenes can be produced by laser- vaporisation of graphite 23|, or en masse by thermal 
methods such as electric-arc discharge 2^ and the combustion of hydrocarbons in a hot 



oxygen-rich environment [25 
ink sticks 



261, 



271] . They have also been found in candle soot [25| , Chinese 
areas of asteroid impact [29| and as byproducts of carbon nanotube synthesis 
301 ]. Despite the technological leaps in the mass production of fuUerenes, the detailed 
atomistic mechanism of fuUerene formation has remained elusive and is hampered by the 
fact that fuUerenes are able to self-assemble out of chaotic environments. Nonetheless, a 
variety of proposed mechanisms exist. Prominent amongst them are: the Pentagon road 
where the fuUerene is created by the addition of C2 units or small carbon particles to the 
dangling bonds of open graphitic cups [31]; the FuUerene road is a similar mechanism to the 
r'entagon road, however the intermediate structures are small closed-cage fuUerenes instead 
32]; Ring Fusion Spiral Zipper where (instead of C2 units) small carbon rings coalesce to 



form a cluster, and then "zip" up and anneal to assemble into a fuUerene |33l . I34| and the 
Shrinking Hot Giant FuUerene road where giant fuUerenes (with attached dangling chains) 
are created via heating a gas of dimers at temperatures 2000-3000 K (the size-up process) in 
a box of side 30 A. These giant fuUerenes are then continuously heated, causing shedding of 
the side chains and evaporation of C2 units which shrinks the giant fuUerenes (the size-down 
process) to sizes as small as the Ceo [35|. Using high-resolution TEM, Huang et al. have 



demonstrated the occurrence of this mechanism where they observed the shrinking of a giant 
fuUerene (C1300) inside the cavity of a multi-walled nanotube. The giant fuUerene shrinks to 
the Ceo, however, below this size, the cage ruptures o pen and evaporates |36|. For detailed 
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reviews of the various formation mechanisms, see 

In this paper, we address the phase transition properties of Ceo and C240 using molecular 
dynamics simulations. We have developed an efficient topologically-constrained pairwise 
forcefield that would allow one to simulate these fragmentation/assembly processes on a 
long time scale (500 ns). The parameters of the forcefield were chosen in order to repro- 
duce the well-known characteristics of the fuUerene molecule. A simple form of the locally 
pairwise interaction allowed us to achieve significant computational gains to observe the co- 
existence of two phases during the phase transition process. The results obtained using the 
developed forcefield are consistent with those obtained with more advanced potentials, see 
Table. IIIII We observe that in a certain temperature range, the fuUerene oscillates between 
two prominently different phases. These states can be considered as solid-like (hollow cage) 
and gas-like phases (gas of dimers) of the fuUerene, and these oscillations correspond to 
the consecutive fragmentation and reassembly of the fuUerene cage at a given temperature. 
We analyze the dependence of the heat capacity of the system on temperature using two 
different approaches: based on a) energy fluctuations in the system and b) differentiation 
of the energy on temperature dependence. Both approaches show that a heat capacity on 
temperature dependence has a prominent peak which is a signature of a first-order-like phase 
transition. 

The paper is structured as follows. In Section [TTl we present the developed forcefield and 
details of the MD simulations conducted. In Section IIIII we present and discuss results of 
the simulations, comparing them to previous works; while in Section IIIID[ we introduce 
entropic corrections to our results by means of a statistical mechanics analysis and show 
the correspondence of this theoretical work to experiments. In Section IIVI we outline the 
conclusions of this work. 



II. THEORETICAL AND COMPUTATIONAL METHODS 
A. Topologically-constrained pairwise forcefield 

The total energy of an A^-particle system is given by the Hamiltonian of the form 



N 



2 



where Fj is the position of atom i, rrii its mass and pj its momentum. The first term 
represents the kinetic energy, while the second is the potential energy based on the forcefield 
V. This forcefield differentiates between two types of interactions in the system: long-range 
van der Waals (vdW) and short-range covalent (cov) bonding. It is expressed as (for all 

N 

V(r,)= Y, yvdw{r,j)+ J2 Vcovirij), (2) 



where rij is the interatomic distance between atoms i and j defined as rij = |rj — r^ |. The sum 
in the covalent term is over the three nearest-neighbors of atom i: atoms k, I and m; while 
the sum in the van der Waals term is over the other (A^ — 4) non-neighbors of atom i. Hence, 
each atom is allocated 3 covalent nearest-neighbors (determined by their positions relative 
to atom i) and A^ — 4 vdW non-neighbors. In this way, the atom is topologically-constrained 
to interact with a certain type of bonding with certain atoms in the system. 

In this work, for computational efficiency we have chosen the Lennard- Jones type of 
potential to treat the interactions between the neighboring atoms although, in principle, 
one could choose different kind of pairwise potentials, such as the Morse potential, to model 
these interactions. 
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where a is the equilibrium distance corresponding to the minimum energy e. For the vdW 
interaction V^dw, the parameters e^dw and cr^dw are based on the nonbonded carbon-carbon 



interaction in 4l|. However, we have adjusted the value of a^dw to 3.0 A from 4.0 A so that 
we could reproduce correctly the bond length of Ceo- 

As for the covalent interaction, there are two types of bonds in the Ceo: the single and 
double carbon bonds. Hence, we have two different equilibrium parameters: 0"^ and cr^ whose 
corresponding minimum energy parameters are e^ and ed- The values of ag and ad reflect the 



equilibrium bond lengths of the single and double bonds in the Ceo |42| , while e^ and e^ are 



)ased on the carbon-carbon bond dissociation energies in ethane and ethylene respectively 



43|. Hence for the Ceo, 



J2 ^^^ ^'^a) = ^^ ('^ik) + Vs (ra) + Vd (ri„), (4) 

j=k,l,m 

since for each atom i, there are three neighbors — two of which form the single bond (also 
known as the 5-6 bond), and the other the double bond (the 6-6 bond). The parameters of 
the forcefield are summarised in Table [B 

TABLE I: The fullerene forcefield parameters 
Type a (A) e (eV) e (kcal/mol) 

Single 1.45 3.81 88.00 

Double 1.38 6.27 144.5 

vdW 3.00 0.00052 0.0121 



In the case of C240, the bonding order is more complicated with bond types intermediate 
between the single and the double bond. To fit the parameters for this (and all other types of 
fullerenes aside from the Ceo), we first calculate all the bond lengths of the DFT-optimized 
C240 and then set these lengths as the equilibrium distance aij for atom i and its neighbor 
j. To calculate the minimum energy parameters, we constructed the following scaling 
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(cd - Cs) + Cs, (5) 



which returns the minimum energy parameters for the Ceo if ctjj = cTs or a^j = a^, corre- 
sponding to (jl]). 

B. Molecular dynamics simulations 

With this forcefield, we have performed constant-temperature molecular dynamics (MD) 



simulations for Ceo and C24o- The time evolution of the system was obtained 



jy integrating 



numerically the Newtonian equations of motion using the leapfrog technique 4J]. The time 



step of the simulation is At=l fs, with each run (at each temperature) being 500 ns long 



(5 X 10^ steps in total). We have assumed the first 500 ps for equilibration purposes, thus 
have only considered data from the remaining trajectory in our analysis. The interaction 
cutoff was set to 15 A and the simulation was performed using periodic boundary conditions 
with a unit cell of length 20 A per side. In contrast to previous simulations that use simple 
10] or the Nose-Hoover thermostat [16], temperature control was achieved 



velocity-scaling 



in this work by using Langevin Dynamics which acts as a stochastic heat bath [45|, |46 |. 
The advantages of Langevin Dynamics are: a) it is known to generate a canonical ensemble 
(unlike velocity-scaling) and b) it simulates the heat transfer of the particles in the heat bath 
through physical "collisions and friction" processes (unlike the Nose-Hoover thermostat) |47J . 
It is described, for each particle i, as 

where 7j is the friction coefficient set to 100 ps~^ and Ri{t) is the random force vector that 
simulates noise in the system. It is defined as a stationary Gaussian process with zero mean: 
{Ri{t)Rj{t')) = 2mi'yikT6{t')6ij, where k is the Boltzmann constant, T the temperature and 
6 is the Dirac delta function. 

The simulation was conducted for temperatures between 3000-8500 K. In the phase 
transition range, we have performed runs in temperature steps of 5 K; while outside this 
range, in steps of 25 K. Simultaneously, we have also varied the e^ parameters of the Cgo 
(es=2.38, 3.25, 3.81, 4.12 and 4.99 eV) to investigate its effect on the phase transition 
temperature. In addition, we have employed a static neighbors list throughout the entire 
simulation, i.e. we have specified the three nearest-neighbors for each carbon atom at the 
beginning of the simulation and it is only with these neighbors that each atom is able to 
'form' and 'break' bonds. We discuss the infiuence of this restriction on the dynamics of the 
system in Sec. IIIIDI 

III. RESULTS AND ANALYSIS 

A. Caloric curves 

We present, in Fig. [1], the temperature dependence of the total energy,as expressed in Eq. 
([T]), for five different e^ parameters (es=2.38, 3.25, 3.81, 4.12 and 4.99 eV). Each of these 
curves can be divided into three distinct parts corresponding to three different phases of 
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FIG. 1: Total energy E of Cgo as a function of temperature for £^=2.38, 3.25, 3.81, 4.12 and 
4.99 eV. Each curve shows a distinct jump in energy corresponding to a phase transition in the 
system. The scattered plot is the time-average total energy and the solid thick line is the cubic 
B-spline interpolation. 

the fuUerene: the region before the abrupt jump in energy (the solid-like phase), the region 
after it (the gas-like phase) and the region of the jump itself (the phase transition). As can 
be observed, an increase in the e^ energy parameter is correlated with an increase in both 
the temperature at which the phase transition occurs and the height of the phase transition 
barrier. This behavior is due to more energy being required to pull apart the bonds in the 
fuUerene. We have also conducted simulations where the e^ energy parameters are varied 
for a given eg parameter which is held constant at 3.81 eV. However, this does not affect 
significantly the thermal behavior of the fuUerene as only 60 out of the 180 bonds in the Ceo 
are double bonds. 

The energy curves corresponding to e^ parameters 3.25 eV and 3.81 eV (equivalent to the 
C-C bond energy in ethylene) show that an energy value of 100-140 eV is required for the 
jump from one phase to another. This result is in very good agreement with that of Horvath 



and Beu who had conducted tight-binding MD simulations and found a multfragmentation 



phase transition occurring for energies between 100-120 eV [17|. It is also in align with 
previous results from other constant temperature |l3l . Il6l | and constant energy |19|, 



48( MD 



simulations which report a result in the range of 90-96 eV. 
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FIG. 2: Snapshots of the fragmentation simulation of Ceo- (a) Deformed Ceo cage at T=4500 K; 
(b) Broken bonds and evaporation of a C2 from the cage at T=5855 K; (c) The broken cage starts 
unravelling and continues evaporating more dimers T=5855 K, 5 ps after the initial C2 evaporation; 
(d) Cage breaks into a web-structure and short chains at T=5855 K; (e) Web-structure decomposes 
to smaller chains and dimers at T=5855 K; (f) Only dimers and short chains at T=6300 K. 



To illustrate the structural changes that occur in the fuUerene, leading to a phase transi- 
tion, we present in Fig. [5] snapshots of the simulation process for es=3.81 eV. The fuUerene 
is found in its solid-Wke phase at T < 5500 K, as demonstrated in Fig. [Uby the linearly in- 
creasing region before the jump. In Fig. [2](a), the fuUerene can be seen to be intact although 
structurally deformed at T=4450 K. At T=5855 K, the carbon bonds are observed to be 
broken and fragmentation occurs with the onset of a C2 evaporation, Fig. El^b). This causes 
the unravelling of the cage structure, as more dimers, short chains and rings are evaporated. 
Fig. [2](c-e). At T > 5855 K, the fuUerene is found to be in its gasASke phase of dimers 



and short chains. As demonstrated in our simulations, the fragmentation/evaporation of 
the fuUerene occurs before the development of a liquid-like 15|, pretzel or linked-chain [l6| 
phases. Instead, the fuUerene is seen to decompose from a stable solid-like cage structure 
to a stable gas-like phase of dimers and short chains, with no stable intermediate structure 
of long carbon chains that would signify a liquid-like phase. This result is also observed by 
Marcos, et al. who attribute the liquid-like phases to insufficient simulation time [l9(]. 



B. The coexistence regime 
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FIG. 3: Top: Bond length as a function of simulation time for a single bond for trajectories with 
es=3.81 eV. At 4450 K (left), the bonds average about 1.5 A. At 5855 K (middle), abrupt jumps in 
the bond length can be seen. These correspond to phase transitions from the solid-like cage phase 
to the gas-like phase of dimers and short chains. At 6800 K (right), the homogeneous distribution 
of bond lengths is characteristic of the gas- like phase. Bottom: Total Energy, E, as a function of 
simulation time over the entire trajectory, including the equilibration process, at T=4450, 5855 
and 6800 K. 



In Fig. [31 we present the bond length of a characteristic single bond in the fuUerene 



as a function of simulation time (for 6^=3.81 eV). The figure spans the whole trajectory 
(500 ns) including the equilibration process which we have taken to be the first 500 ps of 
the simulation. At 4450 K, the system is located in the region before the phase transition, 
Fig. [H and correspondingly, the bond lengths do not fluctuate more than 2.0 A, indicating 
that although the fullerene might be vibrating and deformed in shape, the cage-like structure 
is still present. At 6800 K, the system is located in the region after the phase transition, 
and the homogeneous distribution of the bond lengths is characteristic of a gas-like phase. 
At 5855 K however, one can observe abrupt columns of homogeneous distributions occurring 
for intervals in the simulation. These regions correspond to periods in the trajectory where 
the fullerene has jumped into a gas-like phase. This on-and-off distribution is characteristic 
of a two-phase coexistence regime that has been observed for finite-size clusters (see pLi] and 
references therein). 

In Fig. [31 we plot the corresponding time-dependent total energy for T=5855 K. Here, one 
can observe the abrupt energy jumps between the solid-like phase (~-325 eV) and the gas- 
like phase (~-200 eV) throughout the entire trajectory. This behavior has not been observed 
in previous studies on fullerene fragmentation and thermal stability which we attribute to 
the short time scales of these simulations, as they are usually on the order of 10 ps to 1 ns, 
whereas our simulations are on the order of hundreds of ns. In Fig. HI we present snapshots 
of the transition from a gas-like phase of dimers and short chains, to the solid-like phase of 
a carbon cage. In other words, this is a re-assembly /formation process of the fullerene cage 
where dimers and small chains link up to form a discernible web-like structure that grows 
into the hollow cage. The cage however, is quite deformed as one would expect at such a 
high temperature, nonetheless it is an obvious fullerene cage. 

C. Heat capacity and phase transition temperature 

In Fig. [5l we have plotted the heat capacity of the Ceo obtained from c^ = dE/dT, where 
E is the total energy at temperature T. The solid lines in Fig. [3] were generated by initially 
interpolating the associated energy curve using cubic B-splines (as can be seen in the thick 
solid lines in Fig. [1]) and then by differentiating these interpolated lines. The scattered data 



TABLE II: Entropy change for the different Cg parameters. 



esieV) 



Tpt{K) 



AS 



2.38 
3.25 
3.81 
4.12 
4.99 



3500 
4950 
5855 
6450 
7450 



0.025 
0.024 
0.024 
0.025 
0.027 



points were calculated using the energy fluctuations in the system |49| : 

^ {{E-{E)y) 



(7) 
(8) 



where {E) is the time-average total energy at temperature T and k is the Boltzmann con- 
stant. The maximum of each heat capacity curve denotes the temperature at which the 
phase transition (PT) is defined to occur. We note that in previous papers, the fragmen- 
tation temperature had also been defined differently from this work. They were taken to 



fuUerene cage 
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indicate the onset of bond-breaking 9|, [lOl, lUl] or the dissociation of a C2 unit from the 



48|. These differences are, in part, due to different simulation times 
and the potentials/methods used. In Table IIIIl we have arranged the values of the fragmen- 



tation (or PT) temperature to compare our results with previous works. 

We note that the maximum values of our heat capacity plots are of an order of magnitude 
higher than available fragmentation studies with published heat capacity data |13l . llSl . Il6| . 
This is clearly due to the sharp jump found to occur in the total energy plots due to the long 
simulations times of this work. The almost constant heat capacity distribution across the 
e^ parameters show that the entropy change, AS* from the solid-Yike phase to the ^as-like 
phase is almost constant as 

A5 = 4^, (9) 



To ' 
where AS" is the change in entropy, AE is the change in energy from one phase to the other 

(see Fig. [T]) and Tq is the phase transition temperature. In Table [TTl we have arranged the 





TABLE III: Fragmentation temperature of Ceo 




Model 


T/ 


Ref. 



Zhang, et al. 5500 [10| 

Wang, et al. 5500 [9] 

Xu k Scuseria 5600 [48] 

Kim & Tomanek 4000 [16] 

Kim, Young k Lee 4800 [15] 

Laszlo 6800 [13] 

This model 5855 

entropy change, AS, corresponding to the different single carbon bond energy parameters, 
e^ and the phase transition temperature, Tq. 

D. Statistical mechanics and correspondence to experiment 

The focus of this work is the phase transition dynamics of fuUerene-Uke systems. We have 
constructed a fuUerene model based on a topologically-constrained forcefield which allows 
for gain in the simulation time. To further speed up the simulation process, we have only 
employed a static neighbors list throughout the entire simulation, i.e. we have specified 
the three nearest-neighbors for each carbon atom at the beginning of the simulation and it 
is only with these neighbors that each atom is able to 'form' and 'break' bonds. In other 
words, we have labeled the bonding between atoms throughout the entire simulation. 

While such bond-numbering is reasonable when the system is in the cage structure, it 
is less so when the system moves towards the phase transition region and in its gaseous 
state. This specific bond-numbering prevents the system from forming new bonds. Thus, 
each carbon atom will only seek out its three allocated neighbors that were specified at the 
beginning of the simulation. In other words, our system has been constrained to only one 
particular combination of bonding that would lead to the fuUerene cage. However, for a sys- 
tem consisting of 30 indistinguishable dimers, there should be 30! 2^° possible combinations 
available for the system, all of which would lead to the same cage structure. This correction 
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FIG. 4: Snapshots of the re-assembly process (gas-Uke phase to sohd-hke phase) from the coexis- 
tence behavior of Cgo in the phase transition region at T=5855 K. All times are in ps and denote 
the interval from the first snapshot. 
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FIG. 5: Geo heat capacity plots obtained from c^ = dE/dT for different e^ parameters. The tem- 
perature at the maximum of each heat capacity curve is denoted the phase transition temperature, 
Tpr: 3500 K (2.38 eV), 4950 K (3.25 eV), 5855 K (3.81 eV), 6450 K (4.12 eV), 7450 K (4.99 eV). 
The scatter plots are calculated according to Eq. ([7]) , while the thick solid lines were differentiated 
from the cubic B-spline interpolation from Fig. [TJ 



has to be taken into account in order to specify the correct phase transition temperature. 
Additionally, pressure corrections must also be accounted. In typical electric arc discharge 



experiments, the gas pressure in such setups vary between 13 to 70 kPa [50 



51 



52 



53|; 



however, the pressure in our simulation box is much higher (~0.4 GPa) due to its small 
volume. 

In order to introduce the necessary corrections, let us first analyze the partition of the 
energy in the system at T=5855 K (the phase transition temperature) when it is in the cage 
and gas states (see Table HVT) . The difference in the total potential energy when the cage is at 
T=5855 K and when the system is optimized is equal to 45.8 eV. This value is equal to |NkT 
(where N=60 is the total number of atoms, k is the Boltzmann constant and T=5855 K), 



being in correspondence with the statistical mechanics description of our system. Table IIVI 
shows also that the remaining single bonds energy in the gas phase is equal to -68.82 eV. 
This energy contribution is due to the continuous formation and fragmentation of unstable 
short chains, thus the system in our simulation does not remain entirely in a phase of 30 
dimers. However, in the following statistical mechanics description, we assume that the 
system can reach the gaseous state of 30 entirely unbound C2 dimers at temperatures close 
to the phase transition temperature. 

One can now calculate the entropy change of the system undergoing a phase transition 
with both bond numbering and pressure corrections taken into account. The entropy change 
of the system experiencing a phase transition can be defined as AS = Sgas — Scage, where 
Sgas is the entropy of the gaseous state and Scage is the entropy of the cage structure of 
the fuUerene. Since there are 2^° 30! ways to assemble a fuUerene from 30 dimers, the term 
kln(2^°30!) should be added to the entropy of the cage structure. The entropy of a dimer 
being in the gaseous phase is proportional to the logarithm of the volume accessible to it. 
Therefore the entropy correction to the gaseous state of 30 statistically independent carbon 
dimers in some arbitrary volume (for instance, the volume Vexp corresponding to experi- 



mental conditions) is equal to kin 
30kln 



■Ve, 



\30 



V,i 



or, in terms of concentrations 



30kln 
The similar contribution accounting for the entropy of the cage structure 
within the given volume should be added to Scage- Thus, the total entropy change of the 
system with the corrections has the following form: 



A^ = ASo~k ln(2'''J 30!) + 30A; In 



f^sim 


-fcln 


f^sim 


_ f^exp _ 


_ f^exp _ 



(10) 



-^ is the entropy change of the system ex- 



where k is Boltzmann constant and ASq 

periencing a phase transition from a fullerene cage to the gas state in a simulation box 

of volume 8000 A^. Here the energy variation AE = AEq + PAV represents the energy 



Type 



TABLE IV: Partition of the energy at T=5855 K 



Cage phase energy (eV) 



Gas phase energy (eV) 



Single bonds 
Double bonds 
van der Waals 



-200.18 

-177.51 
+7.045 



-68.82 
-173.67 
+ 2.399 



Total 



-370.65 



-240.09 



difference between the cage structure and the system in the gas phase. We have assumed 
that the short chains of carbon dimers fragment almost at the temperature of the phase 
transition, therefore Aii^o=200.18 eV is the energy difference between the cage state and the 
state where no single bonds are present. This value was obtained from molecular dynamics 
simulations performed at the temperature 5855 K. However, for any arbitrary temperature 
AEq can be calculated with approximately 5% accuracy as follows: AEq(T) = E^h — fkT. 
The first term in the right hand side of the expression Egi, = GOe^ = 228.6 eV is the to- 
tal energy of single bonds for the equihbrated structure (see Table [T]) and the second term 
is the potential part of the vibrational energy of the dimers within the cage. The term 
PAV = {Ngas — Ncage)kTo = 29kTo is the contribution to AE which arises due to the pres- 
ence of pressure in the system. This term accounts for the energy needed to overcome the 
external pressure in the system during the process of fuUerene fragmentation. To =5855 K 
is the phase transition temperature obtained from the simulations while Usim and riexp are 
the concentrations of carbon dimers in the simulation and at experimental conditions re- 
spectively. The second term in the right hand side of Eq. (ITOl) accounts for the permutation 
correction of the entropy in the cage state as discussed above. The third term represents 
the entropy change of the gas of 30 carbon dimers with variation of the system volume from 
Vsim to Vexp, while the fourth term is a similar variation but for the gas of fullerene molecules 
instead. We have assumed that, in the gas phase, the dimers are statistically independent. 
Hence, Eq. (fTOj) can be rewritten as: 



A F 
AS = —^ -k ln(2^° 30!) + 29k In Usim - 29k In 

^0 



P 



where P = n^xpkT is the pressure in the system. Here n~^ = (Vbox — Ka;d)/30, where 
Vfeox=8000 A^ is the volume of the simulation box and Vexd is the excluded volume due to 
the van der Waals repulsion between the atoms. The excluded volume is estimated through 
the distance at which the repulsion energy between dimers is equal to kT. This distance 
r{T) is calculated as a root of the following equation: 



(7) 



12 



(T\' 



.r / 



- kT, (12) 

where e and a are the parameters of the forcefield given in Table [H Knowing the dependence 
of r on temperature and the fact that each pair of atom are not allowed to be closer than the 



distance r, one can calculate the excluded volume: Vexci{T) = 30 ■ ^r'^, which at T=5855 K, 
is equal to Vexci=822 A^. 

Knowing AS as a function of pressure, one can then evaluate the phase transition tem- 
perature as a root of the following expression: 



T 



Esb - IkT ■ 30 + PAV{T) 



(13) 



AS{T) 

where i?s6 = SOe^ = 228.6 eV is the total energy of the single bonds for the equilibrated 
structure (see Table [T]), while the second term in the nominator is the potential part of 
vibrational energy of the dimers in the cage. These two terms describe the energy difference 
between the cage and the gaseous state of the fuUerene at the phase transition temperature. 
The third term PAV{T) = 29kT accounts for the energy necessary for increasing the volume 
of the system during the fragmentation process. The dependence of the phase transition 
temperature on pressure is presented in Fig. [61 
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FIG. 6: Dependance of Cgo phase transition temperature on pressure in the gas phase is plotted 
by solid line. Data points in the figure refer to different carbon-arc experimental conditions: 
Markovic et al. 
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From Fig. [HI it is seen that for the pressure range of 10 100 kPa the phase transition 

temperature is between 3800-4200 K, b eing in a good agreement with the temperatures 



reported in the arc discharge experiments 5l|, ISSj. It is important to note that our statistical 



mechanics model is developed for for an thermodynamically equilibrated system. However, 
in electric arc discharge experiments the condition of local thermal equilibrium for the arc 
region can be justified due to high current density (~320 A/cm^) and high pressures (40- 
70 kPa) |52|. Fig. [U] shows that with lowering the pressure, the phase transition temperature 
decreases. However, at these conditions, the equilibration time increases exponentially with 
the decrease of the phase transition temperature and at some point it should exceed the 
available experimental time limits. 

E. The C240 

In Fig. [7] we have plotted the total energy of C240 as a function of temperature and 
compared it with that for Ceo- As shown, the phase transition of the C240 occurs at a lower 
temperature T ~ 5500 K, while the energy difference between the solid-like and gas-like 
phases of the C240 is ~750 eV — five times larger than that required for phase transition of 
the Ceo- To illustrate this process, we have included the snapshots of the C240 fragmentation 
process in Fig. [HI Similar to the Ceo, the C240 also undergoes deformation of the cage 
which eventually leads to the evaporation of a dimer as the first step towards complete 
disintegration. However, we observe that the C240 is more readily to form larger fragments, 
chains and rings and, unlike the Ceo, it does not reach the coexistence regime at the the phase 
transition temperature. This regime for C240 is expected to happen at a longer simulation 
times. As before, we have derived the heat capacity for the C240, see inset to Fig. [3, and 
from its peak location we have determined the fragmentation temperature being equal to 
be 5500 K. Note that in both figures, the vertical scale on the right-hand side corresponds 
to Ceo while the left-hand side to C240- 

We have demonstrated from our simulations, the fragmentation process of the fullerene 
using our constructed forcefield. Our simulation does not show liquid-like phases of the 
fullerene and the phase transition is seen to occur between the solid-like and (^as-like phases. 
The solid-to-liquid transition are not thermally stable and will fragment into a gas-like phase 

nn 

if the simulation is run long enough p^, 119|]. Most studies in this area have been performed 
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FIG. 7: Total energy of C240 as a function of temperature compared to that for Cqo- Open squares 
correspond to C240 and closed triangles to Cgo- The left vertical scale corresponds to C240 and the 
right to Geo- Inset: heat capacity plots for G240 and Ggo, the horizontal axis is temperature in K 
and the vertical axes are the heat capacity in eV/K. 



with MD simulations that are between 10 ps to 1 ns, although it is known experimentally 
that both laser-induced and thermal fragmentation of fuUerenes occur on the /is timescale 







im. 



24j ]. This therefore implies that since the time scale of any simulation is much 



shorter than the scale of the experiments, the MD fragmentation temperatures will therefore 
overestimate the experimental ones. However, longer simulation time should allow one to 
observe phases that are stable in the fuUerene fragmentation process, and in this field, 
our simulations are two orders of magnitude longer than other MD simulations performed. 
Nonetheless, one must also take into account that in our model, the carbon atoms are 
unable to form bonds with atoms that have not been designated as its neighbors, and this 
might affect both the observation a liquid-like phase and the energetics of the fragmentation 



a) T = 5500 K 
t = 15ns 




b) T = 5570 K 
t = 6.36 ns 




c) T = 5570 K 
t = 60 ns 




d) T = 5750 K 


e) T = 5570 K 


fj T = 7500 K 


t= 121 


ns 


t= 134 ns 


t=15ns 


'^TM 


r 


. "^ '' 




■ 


• 

J 

* 


• 1 . 


• * • V -- ■ ~ ■ 



FIG. 8: Snapshots of the fragmentation simulation of C240: (a) T=5500 K, distorted fullerene cage 
evaporates of a C unit, (b) T=5570 K, broken bonds and formation of hole in cage corresponding 
to further evaporation of C units, (c) T=5570 K, evaporation of a dimer, (d-f) fragmentation of 
the cage follows from the onset of bond-breaking and evaporation until complete dimerization in 
(f) at T=7500 K. Note: All times in the figure correspond to the simulation time at the particular 
temperature only. 



process. However, our result of a phase transition occurring at T=3800-4200 K and requiring 
an energy of 130-140 eV is in good agreement with previous work in this field ly, [l7|, Il9| . 



We have also shown a novel coexistence behavior of the Cgo in the phase transition region 
and have observed the reassembly process of the fullerene cage from a gas of carbon dimers. 

IV. CONCLUSIONS 



We have conducted constant-temperature molecular dynamics simulations of Ceo and 
C240- Without pressure and entropy corrections, our simulations demonstrate that both Ceo 



and C240 undergo a phase transition to a gas-like state at 5855 K and 5500 K respectively. 
For the Cgo we report a two-phase coexistence behavior where the fullerene continuously 
oscillates between a solid-like hollow cage and a ^as-like state of carbon dimers. Such 
oscillation corresponds to continuous fragmentation and reassembly of the fullerene cage (see 
Figs. Eland Hj). We demonstrate that the process of fullerene fragmentation and reassembly 
is a first-order- like phase transition between gaseous and cage states in the finite system. 
The developed topologically-constrained forcefield allowed for simulations that are 500 ns 
long, more than two orders of magnitude longer than those conducted in previous works and 
the results using this forcefield correlate well with those from more sophisticated models 
(See table UTTl) . Using statistical mechanics arguments, we have shown that in the absence of 
bond-numbering and taking into account pressure corrections to correspond to experimental 
conditions, the phase transition temperature of Cgo would be found in the range of 3800- 
4200 K (corresponding to pressures between 10-100 kPa). We have also demonstrated a 
good correspondence of our model to experimental conditions in the arc discharge method 
(See Fig. E]). 
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